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Abstract.  A  new  class  of  schemes  of  the  DSMC  type  for  computing  near-continuum  flows  has  been  recently  suggested:  the 
time-relaxed  Monte  Carlo  (TRMC)  methods.  An  important  step  preceding  the  wide  use  of  these  schemes  is  their  validation  by 
classical  homogeneous  and  one-dimensional  problems  of  gas  dynamics.  For  this  purpose,  a  plane  Couette  flow  is  considered 
in  the  present  paper.  A  comparison  of  TRMC  results  with  the  data  obtained  by  time-proved  schemes  of  the  DSMC  method 
(here  we  used  the  Majorant  Frequency  Scheme)  in  a  wide  range  of  Knudsen  numbers  and  for  different  values  of  wall  velocity 
is  presented. 


INTRODUCTION 

The  Direct  Simulation  Monte  Carlo  (DSMC)  method  [1]  is  the  best  developed  and  powerful  tool  of  simulation  of 
multidimensional  gas  fbws  at  the  kinetic  level.  The  main  drawback  of  this  method  is  the  high  computational  cost 
in  computing  near-continuum  fbws  (with  Knudsen  numbers  Kn  <  0.01).  Therefore,  it  is  of  interest  to  develop  new, 
numerically  effective  schemes  of  the  DSMC  method. 

A  new  class  of  promising  schemes  of  the  DSMC  method  for  computing  near-continuum  fbws  has  been  recently 
suggested  [2]  :  the  time-relaxed  Monte  Carlo  (TRMC)  methods.  An  important  step  preceding  the  wide  use  of  these 
schemes  is  their  validation  by  classical  homogeneous  and  one-dimensional  problems  of  gas  dynamics.  For  this  purpose, 
a  plane  Couette  fbw  is  considered  in  the  present  paper.  The  plane  Couette  fbw  is  a  simple  one-dimensional  fbw  where 
typical  features  of  multidimensional  gas  fbws  are  manifested:  nonuniform  distribution  of  the  gas  in  space,  velocity 
slip  and  temperature  jump  near  the  surface  of  the  plates,  strong  infiience  of  viscosity,  and  heat  transfer  in  the  gas. 

Exact  duplication  of  these  effects  in  the  course  of  simulation  indicates  whether  the  method  can  yield  a  solution 
corresponding  to  the  solution  of  the  Boltzmann  equation.  For  these  reasons,  the  plane  Couette  fbw  is  fairly  suitable 
for  validation  of  new  schemes  of  the  DSMC  method,  such  as  TRMC  schemes. 

The  goal  of  the  present  study  is  to  compare  TRMC  results  with  the  data  obtained  by  time-proved  schemes  (Majorant 
Frequency  Scheme  (MFS))  of  the  DSMC  method  in  a  wide  range  of  Knudsen  numbers  for  different  values  of 
dimensionless  wall  velocity.  In  analyzing  the  modeling  results,  special  attention  will  be  paid  to  high-order  moments  of 
the  distribution  function  (such  as  the  heat-fhx  vector),  which  offer  an  explicit  demonstration  of  the  simulation  quality, 
in  contrast  to  low-order  moments  (such  as  density  and  velocity). 


SPATIALLY  INHOMOGENEOUS  BOLTZMANN  EQUATION 


Let  us  consider  the  Boltzmann  equation,  which  is  ID  in  space  and  3D  in  velocity.  Then  the  equation  reads 


d£ 

dt 


df  1  n(,. 

■”'S  = 


where  the  equation  is  written  in  dimensionless  form  and  Kn  denotes  the  Knudsen  number.  If  the  space  is  divided 
into  a  certain  number  of  cells,  and  the  time  is  divided  into  discrete  intervals  of  size  At,  then  it  is  possible  to  solve 
the  Boltzmann  equation  by  a  splitting  method.  First,  the  collisional  step  is  performed  in  each  cell,  as  in  the  spatially 
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homogeneous  case,  then  the  convection  step  is  computed  for  a  free  fbw  x  ,■  =  x,  +  vv/Af  over  all  particles.  If  the  particle 
crosses  the  solid  surfaces  then  the  particle-surface  interaction  is  performed  (refaction  of  the  particle  with  a  new 
velocity).  The  convection  step  is  performed  in  the  same  manner  for  both  schemes,  whereas  the  implementation  of  the 
collisional  step  is  different  in  MFS  and  TRMC  schemes. 


MAJORANT  FREQUENCY  SCHEMES 


The  DSMC  method  is  traditionally  considered  as  a  method  of  computer  simulation  of  the  behavior  of  a  great  number 
of  model  gas  molecules.  Each  simulated  particle  is  then  regarded  as  representing  an  appropriate  number  Fn  of  real 
molecules.  The  state  of  each  model  particle  is  characterized  by  its  coordinate  r  and  velocity  v.  The  evolution  of 
a  system  of  N  model  particles  can  be  represented  as  a  jump-like  motion  of  a  point  {R,V}  =  {ri,vi, 
in  the  6A'-dimciisional  phase  space.  The  DSMC  method  can  be  then  treated  as  computer  simulation  of  the  6N- 
dimensional  random  jump-like  process.  For  such  a  simulation,  we  represent  the  trajectory  of  the  random  process 
from  the  initial  trajectory  point  of  the  random  process  (to ,  Ro  ,  To),  as  a  transition  in  time  between  subsequent  collisions 
(free  motion)  from  the  state  (fo,RoiFo)  to  the  state  (ti,/?i,Vo),  and  the  changes  in  particle  velocities  after  collisions, 
i.e.,  the  transition  from  the  state  (/q,/?i,Vo)  to  MFS  construction  can  be  illustrated  by  a  simple  case  of  a 

spatially  homogeneous  gas  fbw.  For  a  spatially  homogeneous  rarefied  gas  fbw  case,  there  is  no  dependence  on  spatial 
coordinates  R.  and  the  process  can  be  simulated  only  in  a  3A'-dimensional  space  of  velocities  {V} . 

The  so-called  ’’master”  kinetic  equations  (MKE)  [3],  [4],  [5]  describe  the  behavior  of  an  A'-particlc  gas  model  with 
binary  collisions.  These  linear  MKEs  transform  to  the  nonlinear  Boltzmann  equation  as  N  — »  °°  and  molecular  chaos 
conditions  are  satisfied  (see,  e.g.,  [7]).  It  is  natural  to  use  directly  these  equations  for  constructing  numerical  schemes 
of  the  DSMC  method. 

Using  a  function  w  -  the  probability  density  of  the  transition  of  a  pair  of  particles  from  (vj,v'  )  to  (v;,v/)  -  Kac’s 
MKE  in  a  spatially  homogeneous  rarefied  gas  fbw  case  is  written  as 


^-tfN(t,V)  +  v(V)fN(t,V) 


j]'lJ  -  vt^dv'jdv'j, 


where  the  sum  over  i  <  j  means  summation  over  N (N  —  1  )/2  collision  pairs,  v(V )  is  the  total  collision  frequency, 

V(F)  =  ^'Lat^ij)8ij<°°, 

V  i<j 


g  is  the  relative  collision  velocity,  and  Ot  (gij)  is  the  total  collision  cross  section.  By  introducing  the  majorant  collision 
frequency  [6] 

N(N-l)r  , 

Tm  —  ^  [<? tT/  (g ) ] max  A  V(V) 

,  summation  over  N(N  —  l)/2  collision  pairs  can  be  avoided,  and  MKE  can  be  transformed  to  be  suitable  for 
constructing  a  random  process  that  describes  the  behavior  of  the  A'-particle  gas  model  with  computational  time  linearly 
depending  on  the  number  of  modeling  molecules  N. 

For  a  spatially  homogeneous  case,  the  time  of  the  next  collision  1 1  has  a  probability  density  distribution 
V,„exp{— V,„(fi  —  to)}-  For  the  transition  (ti,V0)  — >  (ti,Vj),  a  collisional  pair  (i,j)  is  uniformly  chosen  from 

N(N  —  l)/2  pairs.  A  real  collision  occurs  with  a  probability  P  =  >  and  with  a  probability  (1  —  P)  the 

velocities  will  not  change,  i.e.,  a  fictitious  collision  occurs. 

For  the  spatially  inhomogeneous  case,  such  a  procedure  is  performed  in  each  collisional  cell,  and  a  detailed 
description  of  MFS  derivation  in  the  case  of  a  spatially  inhomogeneous  gas  fbw  can  be  found  in  [8], 


TIME-RELAXED  MONTE  CARLO  METHODS 

Time-Relaxed  schemes  are  based  on  a  formal  expansion  of  the  solution  of  the  Boltzmann  equation,  called  Wild  sum 
Expansion  ([9,  2]),  and  on  its  subsequent  probabilistic  interpretation.  Note  that  this  probabilistic  interpretation  holds 
uniformly  for  any  time  step. 
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From  the  Wild  sum  expansion  of  the  Boltzmann  equation,  the  following  class  of  numerical  schemes  is  obtained: 


r+l(y)  =  (1  -  T)  £  T*/?(v)  +  T'”+1M(v). 

k= 0 

Here,  f"  «  f(nAt),  At  is  a  small  time  interval,  and  t  =  1  — 

Such  methods  preserve  the  moments  and  are  asymptotically  preserving,  i.e.,  in  the  limit  Kn  — >  0,  the  function  is 
projected  into  a  Maxwellian  function.  The  time  accuracy  depends  on  the  level  of  truncation  m. 

A  generalization  of  the  above-given  schemes  (generalized  TRMC)  can  be  written  as 

m 

r+1(v)=£«(v)+Am+1M(v). 

k= 0 

where  the  weights  A =A^{ t)  are  nonnegative  functions  satisfying  some  conservation  consistency  and  the  asymptotic 
preservation  property.  A  more  detailed  description  can  be  found  in  [2], 

TRMC  schemes  are  a  direct  consequence  of  the  probabilistic  interpretation  of  the  truncated  expansion  of  the  Wild 
sum  expansion.  Here  we  illustrate  the  third-order  TR  method.  Even  if  we  do  not  consider  a  spatially  homogeneous 
problem,  we  can  continue  to  assume  that  /  can  be  normalized  to  unity,  so  that  it  has  the  meaning  of  the  probability 
density. 

We  recall  ([2])  for  the  case  of  the  first-  and  second-order  TR  schemes.  Third-order  TR  scheme  (TRMC3):  For  m  =  3, 
the  generalized  TR  schemes  yield 


fn+1  =A0f  +A\f\  +A2/2  +A3/3  +A4M 

Since  A/,,  k  -  0, 1,2, 3, 4,  are  nonnegative  numbers  that  are  summed  up  to  give  unity,  they  are  probabilities  ([2]). 
Furthermore,  f\  =  P(f"  ,f") /li  is  also  a  probability  density.  The  probabilistic  interpretation  of  the  equation  given 
above  is  as  follows. 

A  particle  extracted  from  f" 

•  does  not  collide  with  a  probability  Ao,  (i.e.,  it  is  sampled  from  /"); 

•  collides  with  another  particle  extracted  from  /"  with  a  probability  A 1  (i.e.,  it  is  sampled  from  the  function  /)); 

•  collides  with  another  particle  extracted  from  f\  (i.e.,  with  a  particle  involved  in  just  one  collision  before)  with  a 
probability  A2  (i.e.,  it  is  sampled  from  the  function  /1); 

•  collides  the  first  time  with  a  particle  sampled  by  /"  and  the  second  time  with  a  particle  extracted  by  f\  (i.e.,  with 
a  particle  already  involved  in  one  collision)  with  a  probability  ^3  (i.e.,  the  particle  is  sampled  from  the  function 

n/1,/1). 

jU  ’ 

•  collides  with  a  particle  extracted  from  /2  (i.e.,  with  a  particle  involved  in  exactly  two  collisions)  with  a  probability 
IA3  (i.e.,  the  particle  is  sampled  from  the  function  F^’/o) . 

j  M 

•  is  replaced  by  a  particle  sampled  from  the  Maxwellian  function  with  a  probability  A4. 

Remark.  In  this  formulation,  the  probabilistic  interpretation  holds  uniformly  in  p At,  in  contrast  to  schemes  based 
on  different  time  discretization  of  the  equation,  which  require  p  At  <  c,  for  some  suitable  c  >  0.  Furthermore,  as 
pAt  — »  °°,  the  distribution  at  the  time  n+  1  is  sampled  from  the  Maxwellian  function.  In  a  spatially  inhomogeneous 
case,  this  would  be  equivalent  to  a  particle  method  for  the  Euler  equations.  In  all  our  numerical  tests,  we  use  TRMC3 
(i.e.  m  =  3). 


NUMERICAL  RESULTS 

The  formulation  of  the  plane  Couette  fbw  problem  is  as  follows.  Let  us  consider  two  infinite  parallel  plates  (L  is  the 
distance  between  the  plates).  The  area  between  the  plates  is  filled  by  a  monatomic  gas  with  an  average  density  n 0. 
The  diffuse  refection  law  is  set  at  the  surface  of  the  plates.  The  plates  have  equal  temperatures  Tw.  The  molecular 
weight  is  m,  the  total  cross  section  is  a,  and  the  molecules  collide  according  to  the  Hard  Sphere  model.  The  left  plate 
is  stationary,  and  the  right  plate  moves  in  the  y  direction  (in  its  own  plane)  with  a  velocity  Vw.  The  mean  free  path  is 


579 


FIGURE  1.  Knudsen  number  Kn  =  0.01,  plate  velocity  S  =  1 


FIGURE  2.  Knudsen  number  Kn  =  1,  plate  velocity  5=1 


Ao  =  1  / (\/2ii()(j).  and  the  Knudsen  number  is  Kn  =  Xq/L.  The  parameter  S  is  the  right  wall  velocity  Vw  normalized  to 
the  most  probable  thermal  velocity  Vt  =  \j2kTw/m  :  S  =  Vw/Vt. 

We  considered  four  cases  summarized  in  the  table  below. 


Case 

Kn 

S 

N  cells 

N  part./cell 

1 

0.01 

1.0 

200 

100 

2 

1.0 

1.0 

20 

100 

3 

0.01 

5.0 

200 

100 

4 

1.0 

5.0 

20 

100 
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FIGURE  3.  Knudsen  number  Kn  =  0.01,  plate  velocity  S  =  5 


FIGURE  4.  Knudsen  number  Kn  =  1,  plate  velocity  S  =  5 


Results  of  the  simulation  are  shown  in  Figs.  1,  2,  3,  and  4.  There  are  profiles  of  density,  temperature,  velocity  of  the 
gas  and  pressure,  shear  stress  tensor  component  Pxx,  and  heat-fhx  component  S  x.  All  the  results  in  the  figures  are 
presented  in  dimensionless  form.  The  density  is  scaled  by  »o,  the  temperature  is  scaled  by  Tw .  and  the  gas  velocity 
is  scaled  by  Vt.  The  pressure  p  and  the  stress  Pxx  are  scaled  by  nomV2,  and  the  heat-fhx  component  S x  is  scaled  by 
HQinV2 .  The  time  step  in  all  cases  was  relatively  small:  At  =  0.1Ax/y/2RT\ . 

The  results  for  Case  1  are  shown  in  Figure  1 .  There  occurs  gas  heating  in  channel,  so  there  is  a  peak  in  the  gas- 
temperature  profile  and  a  minimum  in  the  density  at  the  middle  point  of  the  channel  (X /L  =  0.5).  For  Case  1,  the 
pressure  is  very  close  to  a  constant  value  across  the  channel  (note  that  the  continuum  approach  for  the  Couette  fbw 
provides  a  constant  value  of  pressure).  The  momentum  conservation  equations  imply  that  the  value  of  the  stress 
component  Pxx  is  constant.  Hence,  the  value  of  Pxx  and  its  constancy  across  the  channel  are  the  good  criteria  for 
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validation  of  new  schemes  (see  the  Pxx  profiles).  For  Case  1,  the  heat  fLix  due  to  changes  in  temperature  in  the  channel 
is  quite  small.  It  is  seen  that  the  profiles  of  the  x  component  of  the  heat  fLix  for  the  TRMC  and  MFS  schemes  are  fairly 
close.  Note  that  excellent  coincidence  of  the  simulations  results  for  the  MFS  and  TRMC  schemes  was  obtained  for 
this  case. 

The  results  for  Case  2  are  shown  in  Figure  2.  For  such  a  fbw  (Knudsen  number  Kn  =  1),  the  gas  pressure  is  not 
constant  across  the  channel,  which  is  associated  with  a  nonequilibrium  velocity-distribution  function.  We  can  see  that 
the  value  of  Pxx  is  constant  across  the  channel.  In  Case  2,  we  obtain  quite  good  coincidence  of  the  simulations  results 
for  the  MFS  and  TRMC  schemes. 

In  Figures  3  and  4,  the  results  for  cases  of  a  higher  wall  velocity  ( 5  =  5)  are  presented.  For  Case  3  (Fig.  3,  Kn  =  0.01), 
the  gas  pressure  near  the  walls  (in  the  Knudsen  layers)  becomes  slightly  lower.  Both  schemes  reproduce  this  effect. 
The  mean  pressure  in  the  channel  outside  the  Knudsen  layers  in  the  TRMC  scheme,  however,  is  3%  higher  than  that  in 
the  MFS  scheme.  Both  schemes  reproduce  the  constant  value  of  Pxx.  There  is  a  small  difference  between  the  profiles  of 
the  heat-fhx  component  5  x.  Figure  4  for  Case  4  (Kn  =  1,5  =  5)  shows  that  there  is  good  accordance  between  TRMC 
and  MFS  in  terms  of  density,  temperature,  pressure,  and  heat  fhx  5  x.  There  is  a  little  difference  (less  than  1%)  in  Pxx 
values  for  the  MFS  and  TRMC  schemes.  The  TRMC  scheme  produces  a  small  bowl  in  the  Pxx  profile.  However,  the 
bowl  is  very  small  (about  0.1%),  so  the  results  obtained  by  the  TRMC  scheme  are  acceptable  for  this  case. 


CONCLUSIONS 

The  applicability  of  the  Time-Relaxed  Monte  Carlo  method  for  simulating  rarefied  gas  fbws  is  analyzed  by  com¬ 
parisons  with  the  DSMC  method.  The  comparison  shows  that  the  TRMC  scheme  yields  results  close  to  the  DSMC 
method  for  a  plane  Couette  fbw  for  Knudsen  numbers  1  and  0.01  and  dimensionless  wall  velocity  1  and  5.  The  results 
obtained  show  that  the  density,  temperature,  and  velocity  profiles  coincide  for  both  schemes  in  all  cases.  There  is  a 
noticeable  difference  only  in  the  profiles  of  the  shear  stress  tensor  component  Pxx  for  cases  with  the  higher  value  of 
dimensionless  wall  velocity  5  =  5  (the  difference  does  not  exceed  1%).  A  possible  reason  for  the  difference  in  results 
is  a  high  sensitivity  of  the  TRMC  schemes  to  the  magnitude  of  the  time  step.  This  issue  requires  additional  investi¬ 
gations.  It  is  also  planned  to  perform  a  comparison  between  MFS  and  recently  developed  Recursive-TRMC  schemes 
[10],  which  are  based  on  the  whole  Wild  sum  expansion  of  the  solution,  and  therefore  are  exact  in  time. 
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